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Abstract 

Fuel-cell  dynamics  have  been  investigated  with  a  variable-resistance  board  applied  to  a  high-temperature  polymer  fuel  cell,  operating  on  pure 
oxygen  and  hydrogen  at  atmospheric  pressure  and  about  150  °C.  A  particular  pattern  has  been  identified  in  the  voltage-current  density  diagram.  A 
model  has  been  developed,  using  a  characteristic  equation  to  represent  the  external  circuit  as  a  resistive  load,  possibly  nonlinear  and  time-varying, 
rather  than  assuming  the  current  density  to  be  an  independent  variable.  The  model  successfully  explains  the  shape  of  the  transient  current-voltage 
paths,  including  the  different  time  constants  observed  when  reversing  a  transient. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

The  interest  in  fuel  cells  is  increasing  because  of  the  role  they 
might  play  in  many  power-delivering  utilities  in  various  scales, 
with  environmental  benefits  and  increased  efficiency.  Many  re¬ 
search  papers  are  published,  but  most  of  them  assume  steady- 
state  conditions  for  the  cell. 

Recently,  interest  has  increased  also  for  the  dynamic  be¬ 
haviour  of  fuel  cells,  as  a  topic  in  its  own  merit  and  as  a  pre¬ 
requisite  for  control  analysis  and  controller  design.  However, 
many  models  provided  in  the  literature  [  1-5]  were  not  intended 
for  control  studies,  and  in  that  context  they  are  often  “flawed” 
because  they  consider  current  density  as  an  input  to  the  sys¬ 
tem.  This  is  unrealistic  from  a  process-control  perspective,  be¬ 
cause  in  reality  the  current  is  determined  by  the  characteristics 
of  the  fuel  cell  and  of  its  external  load  [6].  For  example,  it  has 
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been  shown  by  Weydahl  et  al.  [7]  and  Rao  and  Rengaswamy 
[8]  that  the  current  does  not  change  stepwise  during  specific 
transients,  but  instead  follows  a  certain  pattern,  which  may  be 
described  by  a  step  change  followed  by  what  appears  to  be  an 
exponential  relaxation.  A  better  approach  for  models  with  real¬ 
istic  inputs  is  therefore  using  the  load’s  characteristic  as  an  input 
instead. 

Conventional  polymer  electrolyte  membrane  (PEM)  fuel 
cells  employ  an  electrolyte  that  relies  on  liquid  water  to  facilitate 
protonic  conduction,  and  operate  at  a  temperature  ranging  from 
50  to  90  °C,  unless  pressurised.  Careful  control  of  the  heat  and 
mass  balances  is  necessary  in  order  to  avoid  drying  out  the 
electrolyte  membrane  or  flooding  the  gas-diffusion  electrodes. 

Polybenzimidazole  (PBI)  represents  a  new  generation  of 
thermally  resistant  polymer  electrolyte  membranes  with  a  re¬ 
ported  glass  transition  temperature  of  420  °C  [9].  The  conduc¬ 
tivity  of  PBI  in  its  pure  state  is  very  low,  about  10“ 12  Sem-1 
[10,11].  Several  studies  of  the  conductivity  in  acid-doped  PBI- 
membranes  for  fuel  cells  have  been  published  in  the  last  few 
years  [12-16],  The  most  frequently  used  dopant  is  phosphoric 
acid,  as  first  introduced  by  Wainright  et  al.  [17].  Phosphoric 
acid-doped  PBI  fuel  cells  allow  for  an  operating  temperature  as 
high  as  200  °C.  He  et  al.  [  15]  found  the  conductivity  of  PBI  to  be 
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Nomenclature 

a  activity 

A  electrode  area  (m2) 

C  specific  capacitance  (F  m-2) 

D  diffusivity  coefficient  (m2  s-  1 ) 

Erev  reversible  voltage  (V) 

F  Faraday’s  constant  (C  mol-1) 

A g  t  specific  Gibbs  free  energy  of  formation  (J  mol- 1 ) 

;  circuit  current  density  (A  m-2) 

ic  crossover  current  density  (A  m-2) 

z’l  mass-transfer-limiting  current  density  (A  m-2) 

;r  reaction  current  density  (Am-2) 

i'o  exchange  current  density  (Am-2) 

I  current  (A) 

k  reaction  rate  constant  (mol  m-2  s- 1 ) 

L  thickness  of  the  diffusion  layer  (m) 

n  number  of  exchanged  electrons 

N  molar  flux  (mol  m-2  s- 1 ) 

p  pressure  (Pa) 

/•  specific  resistance  (£2  m2) 

R  gas  constant  (J  mol- 1 K- 1 ) 

t  time  (s) 

T  temperature  (K) 

V  voltage  (V) 

x  distance  from  the  reaction  surface  (m) 

Greek  letters 

a  charge-transfer  coefficient 

e  porosity 

p  activation  overvoltage  (V) 

X  equivalent  thickness  of  the  diffusion  layer  (m) 

£  tortuosity 

r  time  constant  (s) 


6.8  x  10-2  S  cm-1  at  200  °C  with  a  FI3PO4  doping  level  of  5.6 
(mole  number  of  H3PO4  per  repeat  unit  of  PBI)  and  5%  relative 
humidity. 

High-temperature  PEM  fuel  cells  based  on  phosphoric  acid- 
doped  PBI-membranes  and  operating  up  to  200  °C  have  several 
advantages  over  conventional  PEM  fuel  cells  working  at  temper¬ 
atures  lower  than  100  °C.  The  increased  temperature  improves 
the  reaction  kinetics  in  the  sluggish  phosphoric  acid  environ¬ 
ment,  simplifies  the  heat  and  mass  transport  management  and 
leads  to  an  increased  tolerance  towards  carbon  monoxide,  which 
is  of  major  importance  for  platinum-based  catalysts.  Li  et  al.  [18] 
demonstrated  the  enhanced  CO  tolerance  in  a  PBI  fuel  cell  at 
elevated  temperatures  and  found  a  CO  tolerance  of  3%  CO  in 
hydrogen  at  current  densities  up  to  0.8  A  cm-2  at  200  °C.  This 
significantly  simplifies  the  pre-treatment  of  the  fuel  and  makes 
the  PBI  fuel  cells  a  viable  alternative  to  the  low-temperature 
PEM  fuel  cells. 

Currently,  there  are  not  many  papers  about  control  of  fuel 
cells.  For  example,  as  of  March  2006,  only  two  articles  about 


fuel  cells  has  been  published  in  the  lournal  of  Process  Control 
[19,20].  However,  sections  about  control  of  fuel  cells  can  be 
found  in  papers  about  dynamics  of  fuel  cells,  or  outside  the 
process-control  literature.  For  instance,  Iqbal  [21]  devised  a 
control  system  for  a  fuel  cell  connected  to  a  wind  power  plant; 
however,  since  he  used  the  Ziegler-Nichols  PID-setting  rules, 
his  resulting  controller  was  too  aggressive,  and  did  not  properly 
reject  disturbances,  especially  the  variation  in  wind  speed.  An 
attempt  to  devise  a  simple  control  system  for  a  PBI  fuel  cell  is 
present  in  a  MSc  thesis  by  Johansen  [22],  but  the  results  were 
not  encouraging;  the  use  of  a  simple  proportional  controller  was 
shown  to  deliver  insufficient  tracking  performance.  Pukrush- 
pan,  in  his  doctoral  thesis  [23],  concentrated  on  the  dynamics 
and  control  of  the  gas  flow  systems,  including  natural-gas 
reformers,  anode  humidifiers,  and  compressors;  however,  based 
on  previous  claims  by  Guzzella  [24],  he  claimed  that  the  time 
constant  of  the  electrochemical  transient  was  in  the  order  of 
magnitude  of  10-9  s,1  which  is  contradicted  by  laboratory 
evidence;  the  electrochemical  transient  is  in  the  range  of 
10-3  to  10  s  for  alkaline  [7,25]  fuel  cells,  and  in  the  range  of 
10-100  s  for  direct-methanol  fuel  cells  [26].  As  it  will  be  shown 
later,  PBI  fuel  cells  have  transients  typically  lasting  between 
10-1  and  1  s. 

Understanding  the  dynamic  behaviour  of  a  PBI  fuel  cell  will 
be  useful  in  developing  efficient  control  algorithms,  which  are 
a  necessary  step  before  applications  using  this  type  of  fuel 
cells  can  be  marketed.  Since  PBI  fuel  cells  require  tempera¬ 
tures  higher  than  the  environment’s,  the  most  likely  applications 
would  be  in  the  higher  range  of  power  ratings,  such  as  the  auto¬ 
motive  and  power-generation  sectors.  The  analysis  presented  in 
this  paper  may  be  adapted  to  other  types  of  fuel  cells,  highlight¬ 
ing  similarities  and  differences  among  different  types  of  fuel 
cells. 

The  objective  of  this  paper  is  to  show  the  development  of 
a  dynamic  model  using  the  external  circuit’s  characteristic  as 
an  input,  and  how  this  model  is  able  to  explain  experimental 
results  observed  on  a  high-temperature  polymer  (PBI-based) 
fuel  cell.  The  first  data  confirming  the  model  came  actually 
from  experiments  already  done  on  alkaline  fuel  cells  [7], 
as  there  are  no  significant  differences  between  this  model 
and  a  model  specifically  designed  for  alkaline  fuel  cells  in 
the  data  range  where  the  experiments  were  carried  out.  The 
measurements  were  later  repeated  on  PBI-based  fuel  cells, 
with  similar  results;  however,  in  the  PBI  fuel  cell,  there  were 
also  superimposing  transients  that  were  found  to  be  much 
slower. 

2.  Experimental  setup  and  procedure 

2.1.  Experimental  setup 

Experiments  were  carried  out  on  a  PBI  membrane-electrode 
assembly,  prepared  in-house  at  Department  of  Materials  Science 


1  In  actuality  Pukrushpan  claimed  10  19  s  in  his  thesis,  but  this  is  most  likely 

due  to  a  typing  error. 
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and  Engineering  of  NTNU,  in  order  to  determine  the  transient 
response  of  the  cell  to  a  variable  load.  The  laboratory  setup  is 
sketched  in  Fig.  1.  This  project  is  part  of  the  FURIM  project 
in  the  European  6th  Framework  Programme,  contract  number 
SES6-CT-2004-502782. 

The  load  itself  was  assembled  as  a  variable-resistance 
board,  with  two  resistances  in  parallel;  in  series  with  one 
of  these  resistances,  a  switch  could  be  opened  and  closed 
manually,  changing  the  total  resistance  of  the  circuit.  The 
value  of  the  resistances  could  be  selected  manually  from  the 
1 1  positions  of  two  corresponding  knobs.  The  resistance  board 
was  recently  developed  by  Helge  Weydahl  in  a  related  study 

[7]. 

The  cell  could  be  kept  at  a  given  temperature  by  a  custom- 
made  external  electric  heater  with  a  feedback  control  loop, 
consisting  of  a  thermocouple  and  an  electronic  control  unit 
(West®  6400  1/16  DIN  profile  controller).  It  should  be  noted 
that,  whereas  the  laboratory  setup  needed  make-up  heat  because 
of  its  small  size,  an  application  will  probably  need  cooling  in¬ 
stead,  since  the  reaction  heat  will  be  sufficient  to  maintain  its 
temperature. 

The  fuel  cell  ran  on  99.5%  hydrogen  and  99.5%  oxygen 
at  atmospheric  pressure.  Oxygen  was  fed  the  cell  directly, 
whereas  hydrogen  was  bubbled  through  a  37.5%  H3PO4 
solution,  in  order  to  make  up  for  any  loss  of  phosphoric  acid 
in  the  cell’s  membrane.  Both  exhaust  streams  were  bubbled 
through  distilled  water,  to  provide  a  simple  visual  flow  feed¬ 
back.  All  streams  were  controlled  with  a  flow-control  unit  by 
Hi-Tek®. 

To  acquire  data,  an  Agilent®  34970A  data-acquisition  unit 
was  used;  the  data  was  later  exported  from  its  interface  program 
into  Octave  to  produce  plots. 

A  polarisation  curve  for  the  cell  was  measured  using  a 
Wenking®  HP  96-20  potentiostat,  whose  voltage  was  set  by  a 
Wenking®  MVS  98  Scan  Generator.  When  measuring  the  polar¬ 
isation  curve,  the  potentiostat  substituted  the  resistance  board  in 
Fig.  1. 


2.2.  Experimental  procedure 

Some  general  conditions  were  maintained  in  all  tests: 

•  the  cell  was  kept  at  a  constant  temperature  of  150  °C; 

•  the  oxygen  flow  was  kept  at  3.7  cm3  s_1; 

•  the  hydrogen  flow  was  kept  at  5.8  cm3  s_1 . 

Since  the  data-acquisition  unit  could  only  make  voltage  mea¬ 
surements,  the  current  through  the  circuit  was  always  calculated 
from  the  voltage  over  a  known,  small  resistance,  integrated  in  the 
resistance  board.  The  voltages  measured  by  the  data-acquisition 
unit  and  the  times  at  which  they  were  taken  were  the  only  mea¬ 
surement  recorded.  In  order  to  have  synchronous  data  for  cur¬ 
rent  and  voltage,  the  time  vectors  associated  with  the  current 
and  voltage  outputs  were  merged;  each  of  the  outputs  was  then 
linearly  interpolated  on  the  merged  time  vector. 

For  variable-resistance  tests,  the  fuel  cell  was  connected  to 
the  resistance  board.  After  having  reached  steady-state  in  cell 
temperature  and  voltage,  the  data-acquisition  unit  was  turned 
on  and  the  board’s  switch  was  flipped.  For  consecutive  tests, 
the  switch  could  be  flipped  again  after  reaching  a  steady-state 
in  voltage.  The  data  could  then  be  saved  as  a  comma-separated 
value  file.  To  obtain  maximum  detail,  the  data-acquisition  unit 
was  set  to  record  data  as  fast  as  possible  (about  every  0.1  s). 

To  obtain  a  polarisation  curve,  the  potentiostat  was  connected 
to  the  fuel  cell  instead  of  the  resistance  board.  Because  of  the 
slow  secondary  dynamics,  the  measurement  took  about  18  h: 
the  potentiostat  cycled  between  0.825  and  0  V  at  a  sweep  rate 
of  25  p. V  s~ 1 .  In  order  to  obtain  a  manageable  number  of  data 
points,  the  data-acquisition  unit  was  set  to  record  data  every  10  s. 

The  gas-diffusion  electrodes  were  prepared  by  spraying  a 
slurry  of  the  relevant  components  with  an  airbrush  (Badger 
No.  100G).  First,  a  layer  of  carbon  support  was  sprayed  onto 
a  carbon  fiber  paper  (Toray,  TGP-H-120),  and  subsequently  a 
catalytic  layer  was  sprayed  on  top  of  the  support  layer.  The 
Pt  catalyst  loading  was  about  1  mg  cm-2  in  total.  The  carbon 


Fig.  1.  The  laboratory  setup  for  the  experiments.  Resistances  Ri  and  R2  can  be  set  to  1 1  different  values  each.  The  fuel  cell’s  voltage  is  measured  by  voltmeter  V\ 
and  passed  to  the  data-acquisition  unit.  Resistance  r,  known  and  fixed,  is  used  to  measure  the  current;  its  voltage  is  measured  by  V2. 
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black  (Vulcan  XC-72)  was  provided  by  Cabot  Corporation.  The 
electrocatalyst,  HiSPEC™  8000  (50%  Pt  on  Vulcan  XC-72)  was 
purchased  from  Johnson  &  Matthey,  while  the  PBI-membrane 
was  provided  by  the  Technical  University  of  Denmark  (DTU) 
as  a  partner  in  the  FURIM  project,  which  is  acknowledged. 
The  membrane-electrode  assembly  (MEA)  used  in  this  work 
was  made  by  hot-pressing  a  sandwich  of  PBI-membrane 
(4cmx4cm)  and  electrodes  (2cmx2cm).  A  more  detailed 
description  of  the  MEA  preparation  is  given  elsewhere  [27]. 

The  MEA  was  tested  in  a  commercially  available  fuel  cell 
(ElectroChem,  Inc.)  with  one  reference  electrode  probe,  dou¬ 
ble  serpentine  flow  field  and  originally  designed  for  low- 
temperature  PEM  fuel  cells. 

To  determine  whether  the  anodic  overpotential  could  be 
neglected,  it  was  measured  by  means  of  a  dynamic  hydrogen 
electrode  in  the  same  way  as  described  by  Kiiver  et  al.  [28]. 
The  dynamic  hydrogen  electrode  current  was  supplied  by  a 
1.6  V  battery,  to  avoid  any  ground  influence.  The  electrode 
potential  was  corrected  for  resistance  by  current  interrupt. 
The  anodic  overpotential  was  measured  with  pure  hydrogen  at 
150°C. 

2.3.  Uncertainty  analysis 

The  preparation  of  the  membrane-electrode  assembly  is 
likely  to  have  a  major  impact  on  the  fuel  cell’s  properties. 
Since  PBI  membranes  are  not  commercially  available  in  stan¬ 
dard  form,  the  variability  encountered  in  reproducing  the  mem¬ 
brane  may  impair  the  reproducibility  of  the  measurements  re¬ 
ported  in  this  paper.  An  effort  was  made  to  make  the  MEAs  as 
reproducible  as  possible,  as  described  in  Seland  et  al.  [27]. 

The  influence  of  oxygen  flow  is  going  to  depend  also  on  the 
number  and  shape  of  the  cathode  channels.  Precision  in  hydro¬ 
gen  flow,  instead,  has  been  found  to  be  of  little  consequence.  The 
performance  of  the  flow  controllers  and  of  the  temperature  con¬ 
troller  will  also  influence  the  resulting  data,  but  this  will  result 
in  systematic  error  rather  than  uncertainty. 

The  precision  with  which  the  board’s  resistances  are  known 
is  going  to  have  a  direct  influence  the  reliability  of  the  measure¬ 
ments  of  current.  Any  unaccounted  nonlinearities  would  also 
disturb  the  measurements. 

Some  transients  had  time  constants  in  the  same  order  of  mag¬ 
nitude  or  even  higher  than  the  sampling  rate.  It  is  possible  that 
some  important  points  might  be  missed  if  a  fast  transient  is  not 
sampled  with  a  sufficiently  fast  sensor,  resulting  in  a  significantly 
different  shape  of  the  transient  measured. 

3.  Modelling  and  analysis 


except  for  diffusion-related  issues,  most  results  will  hold  also 
for  Nafion® -based  fuel  cells  as  well. 

For  the  purposes  of  process  control,  high  modelling  accuracy 
is  usually  not  necessary:  it  is  more  common  to  produce  a  model 
that  captures  the  essential  behaviour  of  the  system,  and  let  a 
feedback  loop  even  out  the  smaller  inconsistencies  as  they  were 
disturbances.  It  is  more  important  to  predict  variations  in  the 
frequency  range  of  interest  than,  for  example,  slow  variations 
that  can  be  taken  care  of  by  a  PI  controller. 

The  model  tries  therefore  to  reach  an  acceptable  compro¬ 
mise  between  the  number  of  parameters  and  accuracy.  Many 
parameters  have  been  lumped  and  some  phenomena  have  been 
neglected  or  simplified,  either  in  order  to  save  computational 
time,  or  because  their  dynamics  were  considered  too  fast  to  be 
of  any  interest  in  a  controller  synthesis. 

3.2.  Assumptions 

The  following  assumptions  have  been  made  in  the  develop¬ 
ment  of  this  model: 

•  diffusion  transients  reach  immediately  their  steady-state,  and 
the  concentration  profile  of  species  is  linear  with  distance 
from  reaction  site  to  bulk; 

•  the  anodic  overvoltage  is  neglected:  Ole  Edvard  Kongstein 
contributed  some  measurements  to  confirm  this  assumption; 

•  transients  are  instantaneous  also  for  the  concentration  of 
species  in  amorphous  phosphoric  acid,  which  depend  on  the 
partial  pressures  of  the  species  at  the  reaction  sites; 

•  the  inverse  (anodic)  reaction  on  the  cathode  is  first-order  with 
respect  to  water  vapour; 

•  the  ohmic  loss  in  the  cell  is  linear  with  current; 

•  cathodic  capacitance  is  constant; 

•  crossover  current  is  constant. 

3.3.  Diffusion 

A  preliminary  analysis  was  run  to  establish  the  influence  of 
diffusion  transients  on  the  response  of  fuel  cells.  Using  actual 
data  and  formulae  as  reported  by  Ceraolo  et  al.  [1],  a  PBI  fuel¬ 
cell  cathode’s  reaction  surface  was  modelled,  with  tortuosity 
|  =  7,  porosity  e  =  40%,  distance  from  reaction  site  to  bulk 
L  =  0.3  mm,  at  a  temperature  of  T  =  150  °C.  These  parameter 
values  are  believed  to  be  severe  enough  to  produce  a  worst-case 
estimate  of  the  time  constant  of  the  diffusion  transient.  The  dif¬ 
fusion  process  involved  water  vapour,  oxygen  and  nitrogen,  and 
was  modelled  with  Stefan-Maxwell’s  equations,  in  the  follow¬ 
ing  form: 


3.1.  Objectives 

The  model  presented  here  was  originally  conceived  as  a  first 
step  in  a  control  analysis  and  synthesis  for  polybenzimidazole- 
based  proton-exchange  membrane  fuel  cells.  This  model  has 
also  been  found  suitable  to  describe  the  transient  behaviour  of 
alkaline  fuel  cells,  when  operating  in  conditions  where  the  dif¬ 
fusion  phenomena  are  not  important  [29];  it  is  expected  that, 


£  dpi 
£2  dx 
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RT 
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( PiNk  -  PkNi) 


(1) 


Since  diffusion  transients  are  caused  by  the  reaction  current 
density  ;r,  they  will  not  in  any  case  be  faster  than  if  s.  If  the 
diffusion  dynamics  is  shown  to  be  faster  than  the  one  associated 
with  ir,  assuming  instantaneous  transients  for  diffusion  is  an 
acceptable  approximation. 
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Furthermore,  the  diffusion  transient  can  have  a  noticeable 
effect  on  the  reaction  rate  only  if  the  partial  pressures  at  the 
reaction  sites  are  sufficiently  low,  that  is,  when  the  cell  is  close 
to  the  mass-transfer  barrier. 

It  will  be  assumed  in  the  following  that  diffusion  transients 
settle  immediately,  and  that  all  partial-pressure  profiles  are  linear 
from  bulk  to  reaction  site.  We  can  then  easily  integrate: 


Pi(x) 


dpt 

Pi.  r  +  -T-* 

OX 


(2) 


Since  the  Stefan-Maxwell  multicomponent  diffusion  Eq.  (1) 
is  valid  for  any  positive  distance  x  <  L  from  the  reaction  site, 
we  can  evaluate  it  with  the  bulk  values  of  partial  pressures, 
Pi ,b  =  pi(L ),  which  are  known;  furthermore,  since  we  neglect 
the  diffusion  transient,  we  can  assume  that  the  flows  of  oxygen, 
nitrogen  and  water  will  be  constant  in  x  and  directly  related  to 
the  reaction  current  ip. 


Nn2  =  0 


(3) 

4  F 

h 

(4) 

IF 

(5) 

Inserting  this  in  the  Stefan-Maxwell  equations  (1)  and  (2), 
the  resulting  direct  formulae  are 


Li;2  RT 

Po2j  =  Po2.b  ~  —Jf 


f  2po2.b  +  pn2o,b 

V  pDo2,n2o 


PNi.b  \  . 

I 
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L%2  RT 

PH2o.T  =  PH2o,b  +  — 


(  2po2,b  +  PH2Q,b 

V  pDo2.H20 
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(7) 


It  can  be  seen  that,  in  this  approximation,  the  reaction-site 
partial  pressures  vary  linearly  with  reaction-current  density.  In¬ 
cidentally,  we  can  obtain  an  approximate  expression  of  the  cur¬ 
rent  b  at  the  mass-transfer  limit,  which  occurs  when  po2.i  0, 
by  rearranging  Eq.  (6): 


II  —  Po2,b 


e  4  F 
L fRT 


(  2po2,b  +  PH20,b 

V  pDo2.n2o 


PN2.b  \  * 
PDo2,N2J 


(8) 


While  this  expression  is  not  an  exact  one,  it  is  very  useful 
as  a  rough,  inexpensive  estimate,  and  as  a  first  guess  in  more 
detailed,  iterative  solution  algorithms. 

Th  2 

InEqs.(6)-(8)thegroup  appears  always  together.  It  is  not 

possible  to  estimate  L,  q  and  s  separately  in  this  model,  so  they 

rt  2 

will  be  integrated  in  a  new  variable,  X  —  X  represents  the 
thickness  of  an  equivalent  diffusion  layer  that  would  result  in  the 
same  reaction-site  partial  pressures  as  in  the  porous  electrode. 


3.4.  Steady-state  activation  overvoltage 


The  model  neglects  any  overvoltages  on  the  anode.  Unless 
differently  specified,  it  will  be  understood  from  now  on  that  all 
electrochemical  terms  refer  to  the  cathode. 


At  steady-state,  the  reaction  current  density  it  is  equal  to  the 
circuit  current  density  i  plus  the  crossover  current  density  ic;  iT 
represents  the  reaction  rate  resulting  from  the  Butler- Volmer 
equation,  i  the  actual  current  flowing  in  the  external  circuit  di¬ 
vided  by  the  cell’s  area,  and  ic  the  internal  losses  such  as  per¬ 
meation  of  reactants  through  the  membrane  or  electronic  con¬ 
ductivity  of  the  membrane. 

Liu  et  al.  [30]  have  recently  studied  the  oxygen  reduction 
reaction  at  the  cathode  of  a  PBI  fuel  cell.  Their  results  (pub¬ 
lished  while  this  paper  was  in  peer  review  and  subsequently  im¬ 
plemented)  indicate  that  the  rate-determining  step  is  first-order 
with  respect  to  oxygen  and  protons: 

02  +  H+  +  e"  -*  02H  (9) 


The  remaining  steps  from  02H  to  H20  are  faster,  and  can 
be  assumed  to  be  instantaneous.  The  oxygen  term  is  the  oxygen 
present  in  amorphous  H3PO4,  and  the  concentration  of  protons 
can  be  promoted  by  water: 

H3PO4  +  H20  ^  H2POJ  +  H30+  (10) 


In  the  data  range  of  interest,  the  data  provided  by  Liu  et  al. 
seems  to  indicate  that  the  concentration  of  H+  is  approximately 
proportional  to  the  partial  pressure  of  water  vapour. 

There  are  no  equivalent  studies  for  the  inverse  reaction  of 
water  oxidation,  so  it  will  be  simply  assumed  to  be  first-order 
with  respect  to  water: 

3IUO  -*  H+  +  e"  +  \02  (11) 


Assuming  that  the  activity  of  oxygen  in  the  amorphous 
H3PO4  phase  is  proportional  to  the  partial  pressure  of  oxygen 
at  the  reaction  site  (with  no  or  negligibly  fast  transient),  the  fol¬ 
lowing  expression  results  for  the  exchange  current  density  [31]: 


i0  =  n  Fk 


a 


(l+cO/2 


(12) 


where  n  is  equal  to  1  and  po  is  the  standard  pressure  (101  325  Pa). 
Finally,  Liu  et  al.  also  provided  an  estimate  for  the  symmetry 
factor  a.  in  the  range  0.42-0.44;  a  —  0.43  will  be  assumed. 

Then,  the  activation  overvoltage  rj  can  be  calculated  for  a 
given  current,  according  to  the  Butler- Volmer  equation: 

ir  =  i0(e“("f/ff7’>''  -  (13) 


Since  this  equation,  which  binds  i]  and  the  reaction  current 
density  iv,  cannot  be  explicitly  inverted  and  expressed  directly  as 
a  function  of  type  rj  =  /(/,-),  an  iterative  algorithm  is  necessary. 

At  this  point,  rj  is  known,  and  there  remains  only  to  calculate 
the  reversible  cell  potential  Elev,  which  can  be  found  from  well- 
known  thermodynamic  data: 

£rev  =  ZML  (14) 

nF 

A g  f  is  calculated  using  the  partial  pressures  at  the  reaction  sites. 
The  ohmic  loss  is  easily  found  as  rceii(i  +  ic).  The  cell  voltage 
is  then  found  as 


V  =  £rev  -  Tj  -  r ceii(i  +  lc)  (15) 

/"cell  is  here  assumed  to  be  constant. 
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3.5.  Dynamics 

During  transients  the  reaction  current  density  z'r  is  different 
from  the  sum  of  current  density  (due  to  the  external  load)  and 
crossover  current  density  (due  to  losses).  Looking  at  the  follow¬ 
ing  differential  Eq.  (1),  it  is  clear  that  the  difference  i  +  ic  —  iv 
is  actually  the  driving  force  of  the  transient  behaviour  of  the 
activation  overvoltage: 


Table  1 


Summary  of  the  main  equations  used  in  the  model 


Overvoltage  differential 

.  _  i+ic-ir 

"  c 

equation 

Butler- Volmer  equation 

ir  =  i0  (yJnF/RT)r,  _  e-(l -a)(nF/RT)^ 

Exchange  current  density 

i0  =  nFk(P^Y 

\  Po  )  \  Po  J 

(12’) 

Oxygen  diffusion 

P02,r  =  /Or,  Pi, b,  A  T) 

(6’) 

Water  diffusion 

PH20,r  =  gOr,  Pi, b,  V  T) 

(7’) 

This  is  the  only  differential  equation  that  we  consider  in  order 
to  compute  the  transient  of  the  fuel  cell.  As  already  mentioned, 
all  diffusion  transients  will  be  neglected;  any  capacitance  or 
inductance  in  any  part  of  the  load  is  also  neglected.  Capaci¬ 
tance  C  is  assumed  to  be  constant;  however,  in  real-world  oper¬ 
ation,  C  varies:  in  PEM  fuel  cells,  a  variation  of  capacitance  as 
a  function  of  cell  potential  has  been  investigated,  for  example, 
by  Parthasarathy  et  al.  [32]. 

The  load  in  the  model  is  generally  allowed  to  vary  with  time, 
and  is  described  by  a  function  of  the  form  I  =  /(V,  t).  It  has 
been  chosen  to  use  this  form  instead  of  V  —  g{I,  t )  because  cur¬ 
rent  I  through  MOSFETs  flattens  out  beyond  a  certain  value  of 
the  source-drain  voltage  V  (depending  on  the  MOSFETs  gate 
voltage),  and  this  yields  an  infinite  number  of  values  of  V  for  a 
single  value  of  I.  MOSFETs  are  interesting  candidates  as  con¬ 
trol  actuators  in  linear  control,  and  therefore  it  is  considered 
important  that  the  model  be  able  to  simulate  them. 

The  differential  Eq.  (16)  describing  the  activation  overvoltage 
looks  simple:  i  is  the  actual  current  in  the  circuit,  whereas  iT 
is  the  reaction  current,  or  the  current  that  corresponds  to  the 
consumption  of  reactants  in  the  fuel  cell.  The  crossover  current 
ic  represents  current  that  is  generated,  but  is  lost  inside  the  fuel 
cell.  The  net  sum  of  these  quantities  is  the  rate  of  accumulation 
of  charge  in  the  model’s  capacitor,  which  represents  the  sum  of 
an  electrostatic  double  layer  capacitance,  a  pseudocapacitance 
due  to  charge  stored  by  adsorbed  and  intermediate  products  at 
the  surface,  and  the  charge  stored  in  the  diffusion  layers  [33]. 
Divided  by  the  capacity  C,  this  sum  yields  the  time  derivative  of 
the  cathodic  activation  overvoltage,  f). 

The  model  can  be  expressed  in  a  diagram,  as  shown  in  Fig. 
2.  There,  one  can  see  the  internal  resistance  rcen,  the  reversible 
potential  Elev,  the  capacitor  C  and  its  overvoltage  rj.  The  bipole 
in  parallel  with  the  capacitor  is  a  nonlinear  voltage-controlled 


current  generator,  which  enforces  the  Butler- Volmer  law  impos¬ 
ing  a  current  ir  as  a  function  of  ?;  (and  many  other  parameters, 
as  per  Eq.  (13)).  This  bipole  is  often  also  modelled  as  a  nonlin¬ 
ear  resistance  [34],  or  a  linearised  value  in  an  area  of  interest. 
Incidentally,  the  ideal  generator  Eiev  and  the  resistance  rceii  in 
series  with  it  can  be  viewed  as  a  Thevenin  equivalent  circuit. 

The  actual  integration  of  Eq.  (16)  is  not  very  simple,  since  zr 
is  a  function  of  many  parameters,  including  rj  itself,  through  the 
Butler- Volmer  equation  (13).  Furthermore,  the  exchange  cur¬ 
rent  density  i o  in  that  equation  depends  on  the  reaction-site  con¬ 
centrations  of  the  reactants  (Eq.  (12)),  which,  in  turn,  depend 
on  their  consumption  rate,  proportional  to  zr.  An  iterative  loop, 
consisting  of  Eqs.  (6),  (7),  (12)  and  (13),  is  therefore  necessary 
to  calculate  zr  at  all  integration  steps,  and  this  is  the  major  com¬ 
putational  cost  of  the  simulation.  This  phenomenon  could  be 
approximated  assuming  a  constant  f'o,  but  a  separate,  empirical 
term  for  mass-transport  overvoltage  should  then  be  introduced  in 
Eq.  (15).  The  complete  set  of  equations  to  solve  is  summarized 
in  Table  1 . 

The  overall  function  that  calculates  the  transient  simulation, 
based  on  the  selected  time  steps  and  the  external  load’s  function, 
calculates  first  the  initial  state  of  the  fuel  cell,  which  is  repre¬ 
sented  by  its  activation  overvoltage  ijq.  To  find  it,  the  model 
assumes  that  the  system  is  initially  at  steady-state. 

3.6.  Time  constants 

Given  the  state  x  of  a  dynamic  system,  a  time  constant  is  the 
factor  r  in  the  following  differential  equation: 

r.t  =  —  x  (17) 

In  this  equation  forcing  terms  were  neglected  for  sake  of  sim¬ 
plicity.  We  will  now  show  how  to  obtain  a  simple  relationship  for 
the  time  constant  of  the  electrochemical  transient  of  fuel  cells. 

The  external  characteristic  is  assumed  to  be  a  generally  non¬ 
linear,  purely  resistive  load  of  the  type  V  —  fifiz).  The  reversible 
voltage  £’rcv,  the  cell’s  internal  resistance  rceii  and  the  crossover 
current  ic  are  assumed  to  be  constant.  The  time  derivative  of  the 
cell’s  voltage  is  then  expressed  as 

dV  d  „  d 

~T~  =  -t\E  ~  *1  ~  r cell  O'  +  *c)]  =  --7-07  +  r cellO  (18) 

at  at  at 

Since  the  operating  point  has  to  be  on  the  load’s  characteristic 
as  well  as  on  the  cell’s  instantaneous  one,  it  can  be  written  that 

dV  d'F  dfl>  dz 

—  =  —  -  —  —  (19) 

dr  dr  dz  (  dr 
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Inserting  in  Eq.  (18)  and  rearranging 


dr 

d7  1  ^eH  + 


''cell  T~ 


dflr 

\  dr] 

di 

i)  dr 

dflr 

y' i+ 

di 

c 


The  overvoltage’s  time  derivative  is 


dr? 

dr 


i  +  rc 


C 

Rearranging: 


di) 
d  zr 


dir 

~dt 


-I 


i  +  Z'c 


C 


Subtracting  Eq.  (22)  from  Eq.  (20),  we  obtain  the  following 
relationship: 


Fig.  3.  Experimental  polarisation  curve  of  a  PBI  fuel  cell  at  150  °C,  measured 
with  a  sweep  rate  of  25  |xV  s_  1 .  The  arrows  indicate  the  direction  of  the  hysteresis 
cycle. 


1/T 

X  (i  +  i'c  -  i'r)  (23) 


However,  given  how  slowly  the  polarisation  curve  was  sampled, 
this  is  unlikely  to  be  the  only  reason.  Likewise,  temperature 
transients,  caused  by  the  varying  power  output,  should  not  have 
influenced  such  a  slow  sampling.  Other  possible  causes  include 
changes  in  the  access  of  oxygen  to  the  reaction  site  or  changes 
in  the  catalyst  during  the  cycle. 


Assuming  that  C  and  rceu  are  constant,  that  functions  'F(z)  and 
z?(zr, . . .)  are  given,  the  time  “constant”  r  depends  essentially  on 
i  and  z'r: 


x(i,  ir)  = 


C 


( 


f cell  + 


(24) 


The  reasons  for  using  i  +  ic  —  zr,  which  represents  the  current 
inflow  in  the  capacitor  and  is  therefore  proportional  to  the  time 
derivative  of  the  overvoltage  r?,  as  a  state  of  the  system,  instead  of 
the  more  natural  r]  itself,  are  that  it  is  possible  to  find  expression 
(23)  explicitly,  without  solving  numerically  the  Butler- Volmer 
equation  (13),  and  that  it  is  not  necessary  to  add  a  forcing  term 
to  account  for  external  effects,  since  they  are  already  included 
through  VP;  furthermore,  the  term  i  +  ic  —  ir  always  returns  to 
0  at  the  end  of  a  transient,  since  it  is  the  transient’s  driving 
force. 


4.  Results  and  discussion 

4.1.  Polarisation  curve 

The  polarisation  curve  obtained  at  the  conditions  described 
in  the  previous  section  is  shown  in  Fig.  3.  It  can  be  seen  that  the 
polarisation  curve  shows  a  fairly  large  hysteresis.  This  behaviour 
has  been  observed  previously  in  other  experiments  on  similar 
PBI  cells  in  Liu  et  al.  [30]. 

Since  Liu  et  al.  showed  that  the  cathodic  reaction  is  first-order 
with  respect  to  the  concentration  of  H+  ions,  which  depends 
on  the  partial  pressure  of  water  vapour,  it  might  be  tempting 
to  blame  the  hysteresis  on  the  lag  in  the  formation  of  water. 


4.2.  Anodic  overvoltage 

The  anodic  overpotential  at  1  A  cm-2  was  measured  to  be 
40  mV,  indicating  a  small  anodic  overpotential  compared  to  the 
cathodic  one.  This  justifies  the  assumption  of  neglecting  the 
anodic  overpotential. 

4.3.  Resistance  steps 

The  results  obtained  by  switching  from  an  initial  resistance 
to  a  smaller  one  and  back  are  shown  in  Figs.  4  and  5.  It  can 


Current,  A 


Fig.  4.  Experimental  operating  points’  path  in  the  V-/  diagram  after  step  changes 
in  the  external  circuit’s  resistance,  starting  from  82.25  £2  and  switching  back. 
The  arrows  indicate  the  path  in  the  I-V  plan. 
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Fig.  5.  Zoom  on  an  experimental  transient  path  for  a  step  change  in  the  circuit 
resistance. 

be  seen  how  the  operating  point  moves  from  one  resistance’s 
characteristic  to  the  other  along  almost  parallel,  straight  lines. 
The  reason  why  the  lines  are  not  exactly  parallel  might  be  that 
sampling  is  not  continuous,  and  the  operating  point  might  have 
moved  from  the  initial  point  on  the  new  characteristic  by  the 
time  a  new  sample  is  measured. 

As  a  reference,  the  steady-state  characteristic  is  plotted  as 
well.  The  last  transient  is  not  exactly  on  the  expected  character¬ 
istic  of  its  external  resistance,  most  likely  due  to  imprecision  in 
the  calculation  of  the  resistance  value. 

It  is  clear  that  some  transients  are  faster  than  others  in  the  de¬ 
scent  to  steady-state;  this  can  be  seen  by  the  number  of  markers, 
roughly  spaced  by  0.06  s,  on  the  trajectories  in  the  V—I  plane.  It 
appears  that  the  electrochemical  transient  was  faster  the  lower 
the  resistance. 

Zooming  in  on  a  specific  transient  provides  further  informa¬ 
tion.  In  Fig.  5  the  second  transient,  the  one  to  5.47  £2,  is  rep¬ 
resented  in  the  V-I  plane;  the  marks  on  the  line  represent  the 
sampling  points.  It  can  easily  be  seen  that  the  descent  to  steady- 
state  on  the  characteristic  of  the  5.47  £2  resistance  was  faster  than 
the  ascent  to  steady-state  on  the  82.25  £2  characteristic  when  the 
switch  on  the  resistance  board  was  turned  off. 

There  is  an  indication,  in  other  words,  that  the  electrochem¬ 
ical  transient  can  exhibit  different  time  constants  depending  on 
both  the  start  and  the  end  point.  This  behaviour  could  be  related 
to  the  phenomenon  observed  by  Johansen  [22],  where  the  same 
proportional  controller  performed  too  aggressively  in  some  op¬ 
erating  ranges,  and  too  mildly  in  others. 

The  values  assumed  by  the  voltage  are  plotted  versus  time  in 
Fig.  6.  The  experiment  is  the  same  as  in  Fig.  4.  It  can  be  seen 
that,  after  a  transient  in  the  order  of  magnitude  of  1  s,  a  slower 
transient  of  generally  smaller  amplitude  appears. 

Wang  and  Wang  [35  ]  measured  a  similar  behaviour  for  a  low- 
temperature  PEM  fuel  cell,  which  had  also  been  described  by 
Ceraolo  et  al.  [  1  ].  In  both  cases,  the  authors  noted  that  slow  tran¬ 
sients  in  proton  concentration  are  the  cause  of  these  transients. 
Since  protons  in  the  amorphous  H3PO4  phase  are  important  in 
PBI  fuel  cells  too  [30],  the  mechanism  could  be  similar  in  prin¬ 


Fig.  6.  Experimental  voltage  transients  after  the  step  changes  in  the  external 
circuit’s  resistance,  starting  from  82.25  £2  and  switching  back. 

ciple.  Compared  to  Wang  and  Wang’s  measurements,  transients 
in  PBI  fuel  cells  appear  to  be  slower. 

It  can  be  seen  from  Figs.  4  and  5  that  all  the  transients  to 
lower  values  of  resistance  seem  to  reach  first  the  lower  branch 
of  the  polarisation  curve,  to  rise  again  to  the  upper  one  as  pro¬ 
ton  concentration  increases.  The  transients  back  to  the  larger 
resistance,  however,  are  much  larger  than  the  gap  between  the 
branches  of  the  polarisation  curve;  they  have  larger  overshoots 
for  transients  arriving  from  higher  values  of  current,  possibly 
because  of  higher  concentration  of  protons.  This  would  seem  to 
strengthen  the  case  for  a  role  of  water  in  the  hysteresis  of  the  po¬ 
larisation  curve  of  Fig.  3,  but,  on  the  other  hand,  these  transients 
are  much  faster  (in  the  order  of  10  s)  than  the  sampling  rate  of 
the  polarisation  curve,  which  would  have  had  enough  time  to 
reach  the  steady-state. 

These  slower  transients  are  not  considered  to  be  an  important 
issue,  since  in  a  control-oriented  perspective,  these  secondary 
transients  can  easily  be  controlled  by  a  simple  PI  feedback  con¬ 
troller  [36].  Therefore,  we  shall  turn  our  attention  to  the  first 
moments  of  the  transient. 

4.4.  Effect  of  flow  disturbances 

It  was  noted  that  just  changing  the  oxygen  flow  rate  would 
change  the  open-circuit  voltage  significantly,  as  shown  in  Fig.  7. 
It  was  initially  supposed  that  variations  in  temperature  could  be 
the  cause  of  this  behaviour:  in  the  experimental  setup  sketched  in 
Fig.  1 ,  gases  enter  the  cell  at  ambient  temperature,  and  could  have 
cooled  the  cathode  to  a  temperature  lower  than  that  measured 
by  the  thermocouple.  However,  subsequent  experiments  with 
gas  pre-heating  showed  that  the  difference  was  minimal. 

The  anode  turned  out  to  be  less  sensitive  to  inlet  hydrogen 
flow.  In  fact,  changing  the  hydrogen  flow  would  not  cause  any 
measurable  disturbance  to  the  operating  point. 

It  can  be  supposed  that  variations  in  proton  concentration  at 
the  cathode  are  the  cause  of  these  transients  too.  When  oxygen 
flow  is  increased,  there  is  an  immediate  increase  in  voltage  due 
to  higher  oxygen  partial  pressure,  and  a  subsequent  decrease 
due  to  reduction  in  proton  concentration,  caused  by  the  reduced 
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Fig.  7.  Experimental  measure  of  the  effect  of  a  step  in  oxygen  flow  from  3.7 
to  18  cm3  s-1  at  150  °C  on  open-circuit  voltage.  The  first  step  is  the  moment 
the  oxygen  flow  is  started  at  the  beginning  of  the  experiment.  The  second  step, 
at  about  t  =  2200  s,  is  the  moment  oxygen  flow  is  stepped  up;  the  last  step,  at 
about  t  =  3600  s,  is  when  oxygen  flow  is  stepped  back  down. 

water  partial  pressure  at  the  reaction  site.2  Conversely,  when 
oxygen  flow  is  reduced,  there  is  a  sharp  reduction  in  voltage  due 
to  lower  oxygen  partial  pressure,  followed  by  a  (this  time  faster) 
recovery  of  proton  concentration. 

This  asymmetry  in  the  dynamics  of  proton  concentration  has 
already  been  noted  in  Ceraolo  et  al.  [1], 

5.  Modelling  results  and  discussion 

5.1.  Diffusion  modelling 

The  result  of  diffusion  modelling  is  shown  in  Fig.  8.  This 
model  predicts  transients,  without  overshoots,  in  the  order  of 
10-1  s  in  case  of  steps  in  the  reaction  rate.  Such  steps  are  im¬ 
possible  in  reality,  since  iT  varies  continuously  as  a  function  of 
As  long  as  the  time  constants  of  the  transients  of  ir  are  larger  than 
10“ 1  s,  the  diffusion  transients  may  safely  be  assumed  instan¬ 
taneous;  they  may  however  be  more  important  where  transients 
in  ir  are  faster. 

Since  we  have  seen  that  transients  usually  take  longer  than 
this  to  complete,  it  is  safe  to  neglect  diffusion  dynamics.  How¬ 
ever,  this  is  only  valid  as  long  as  transients  in  overvoltage  are 
not  faster  than  the  ones  in  diffusion. 

5.2.  Model  performance 

The  model,  implemented  in  C++,  can  calculate  a  resistance- 
step  transient  in  a  few  seconds  of  computation  time  on  a  2.4  GHz 
Pentium4  computer  under  Linux.  Previous  experience  with  lan¬ 
guages  such  as  Scilab  required  about  1  min  for  the  same  tran¬ 
sients,  with  more  marked  difficulties  when  integrating  close  to 
the  mass-transfer  limit. 


2  Even  if  we  are  at  open  circuit,  crossover  current  is  causing  some  reaction  to 
take  place,  and  some  water  to  be  formed. 


Fig.  8.  Simulated  concentration  transient  on  the  cathode  of  a  PBI  fuel  cell;  the 
current  steps  from  0  to  800  A  m  2  at  t  =  0.  The  upper  profile  is  oxygen,  the 
lower  water  vapour. 


Current  density,  A  cm'2 

Fig.  9.  Plot  of  residuals  (sample  minus  model  values)  for  the  model  applied  to 
the  upper  branch  of  the  polarisation  curve  shown  in  Fig.  3. 


Current,  A 

Fig.  10.  Simulated  transient  trying  to  replicate  the  experiment  in  Fig.  5.  The 
value  assumed  for  the  capacitance  is  1 500  Fm~2 .  The  markers  are  spaced  by 
0.06  s,  except  those  at  switching  points  that  are  spaced  1  p,s. 
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Table  2 


The  results  of  the  parameter  regression  for  parameters  k,ic,  A  and  rcen  with  the  data  of  the  upper  branch  of  the  polarisation  curve  presented  in  Fig.  3 


Estimate 

d=cr 

Units 

Correlation  matrix 

k 

115.7 

55.9 

(xmol  m-2  s-1 

i 

-0.590 

0.672 

0.847 

l 

25.82 

0.349 

mm 

-0.590 

1 

-0.740 

-0.648 

ic 

33.63 

93.94 

Am-2 

0.672 

-0.740 

1 

0.518 

r cell 

26.07 

5.96 

(jl  £2  m2 

0.847 

-0.648 

0.518 

1 

5.3.  Parameter  regression 

A  parameter  regression  has  been  carried  out  to  verify  whether 
and  how  well  the  proposed  model  fits  the  results.  We  were  ini¬ 
tially  sceptical  about  a  parameter  regression  to  find  k  and  other 
parameters,  because  of  the  large  hysteresis  of  the  polarisation 
curve  in  Fig.  3:  however,  we  noticed  the  model  fits  very  well  the 
upper  branch  of  the  curve,  and  were  able  to  regress  to  some  satis¬ 
factory  parameter  values  and  produce  a  residual  plot.  However, 
the  lower  branch  does  not  qualitatively  fit  even  with  different 
values. 

The  residuals  are  plotted  in  Fig.  9,  where  it  can  be  seen  that 
the  model  predicts  the  cell’s  voltage  in  the  upper  branch  with 
an  approximation  of  ±  10  mV  until  it  reaches  the  mass-transport 
limit.  The  results  of  the  parameter  regression  are  shown  in  Table 
2.  Estimates,  standard  deviations,  units  and  the  parameters’  cor¬ 
relation  matrix  are  given. 

There  is  significant  correlation  between  internal  resistance 
rcen  and  the  reaction  constant  k,  since  in  most  experimental 
points  the  effect  of  increased  resistance  is  similar  to  that  of  a 
reduced  reaction  constant.  Similarly,  there  is  high  correlation 
between  crossover  current  density  ic  and  the  effective  thickness 
of  the  diffusion  layer  X:  since  information  on  X  is  gathered 
mostly  at  data  points  close  to  the  mass-transport  barrier,  if  ic 
is  larger,  X  must  be  thinner,  to  allow  more  reaction  current 
(;r  =  i  +  ic)  to  pass  at  the  mass-transport  limit  in  order  to 
result  in  the  same  measured  values  of  i.  To  reduce  these 
correlations,  one  may  gather  data  points  from  conditions 
where  the  effects  of  these  parameters  are  going  to  vary 
independently,  such  as  different  values  of  temperature  or  gas 
composition. 

The  model  does  not  fit  the  lower  branch  of  the  polarisation 
curve.  It  is  possible  that  the  lower  branch  follows  different 
kinetics,  or  other  factors,  such  as  catalyst  deactivation  and 
transients  in  proton  concentration,  interfered  in  the  measure¬ 
ment  and  caused  the  hysteresis  observable  in  the  slow  scan 
(Fig.  3). 

Some  computational  difficulties  were  experienced  as  the  re¬ 
gression  algorithm  could  sometimes  converge  to  a  local  mini¬ 
mum  different  from  the  global  one. 

It  was  not  possible  to  perform  a  proper  parameter  esti¬ 
mation  for  the  cathodic  capacitance,  as  the  slower  secondary 
transients  shown  in  Fig.  6  are  not  well  understood,  and  would 
have  interfered  in  the  regression  procedure.  Furthermore, 
the  exact  time  of  switching  could  not  be  logged  in  the 
measurements.  However,  a  value  in  the  range  of  l-2kFm~2 
seems  to  reproduce  the  observed  behaviour  with  acceptable 
approximation.  The  value  1500  Fm-2  has  been  used  in 
Fig.  10. 


5.4.  Path  on  the  V -I  diagram 

The  cause  of  the  sudden  iso-??  “kick”  in  the  dynamic  path 
diagram,  as  can  be  seen  in  the  simulation  in  Fig.  10  and  in  ex¬ 
perimental  results  in  Figs.  4  and  5,  is  caused  by  the  simple  fact 
that  the  operating  point  must  always  lie  on  the  load’s  character¬ 
istic,  when  the  load  is  purely  resistive.  It  is  possible  to  determine 
to  which  coordinates  on  the  load’s  characteristic  the  operating 
point  will  jump,  considering  that  the  activation  overvoltage  ?? 
varies  continuously  with  time;  when  a  step  change  in  the  outer 
load  happens  at  time  to,  r)  will  be  the  same  for  times  f("j~  and  t^ . 

For  every  point  in  the  steady-state  polarisation  curve  we  can 
identify  its  iso-  ??  line,  which  represents  the  set  of  operating  points 
to  which  the  fuel  cell  can  jump  instantaneously.  This  line  is  the 
characteristic  of  the  cell’s  internal  resistance.  The  equation  is 
the  same  as  Eq.  (15),  but  the  first  two  terms  (£'rev  and  ??)  are 
now  constant  for  each  instantaneous  characteristic,  since  they 
depend  on  iT. 

V(I)  =  Eiev  -  n  -  reel,  f f  +  ic)  (25) 

constant 

We  can  visualise  this  line  as  a  sort  of  time-varying  char¬ 
acteristic :  when  the  external  load  changes,  the  operating  point 
moves  immediately  to  the  new  intersection  between  the  new 
load’s  characteristic  and  the  iso-??  line  (the  instantaneous  char¬ 
acteristic).  As  the  value  of  i  will  change  accordingly,  i  +  ic  —  ?r 
will  no  longer  be  equal  to  zero,  and  this  will  cause  the  overvolt¬ 
age  rj  to  vary  with  time  according  to  differential  Eq.  (16).  This 
means  that  the  iso-??  line  will  start  moving  vertically  according  to 
the  characteristic  Eq.  (25),  until  the  intersection  between  it  and 
the  load’s  characteristic  is  again  on  the  steady-state  polarisation 
curve,  where  by  definition  i  +  ic  =  iT. 

Using  the  load  characteristic  as  the  model’s  input  has  allowed 
to  simulate  real  laboratory  results,  and  has  given  the  possibility 
of  simulating  the  implementation  of  MOSFETs  as  manipulated 
variables  in  control  of  fuel  cells.  MOSFET  transistors  are  loads 
whose  characteristic  can  be  changed  by  manipulating  a  gate 
voltage.  The  load’s  current  enters  from  a  source,  and  exits  from 
a  drain;  the  gate  is  a  third  connection  point,  whose  voltage  com¬ 
pared  to  the  source  is  the  manipulated  variable;  since  no  net 
current  flows  from  the  gate  to  either  source  or  drain,  little  or  no 
power  is  consumed  by  the  process  of  changing  the  MOSFET’s 
characteristic.  MOSFETs  can  be  imagined  as  fail-closed  fluid- 
flow  valves:  when  no  gate  voltage  is  applied,  no  current  may 
pass.  As  the  gate  voltage  increases  beyond  a  certain  threshold, 
the  characteristic  detaches  from  the  open-circuit’s  and  allows 
more  current  to  pass.  Normally,  for  every  gate  voltage  there  is  a 
maximum  current  that  can  pass  through  the  MOSFET,  regard¬ 
less  of  the  voltage  across  the  source  and  drain  of  the  transistor. 
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More  details  on  transistors  and  MOSFETs  in  particular  can  be 
found  in  textbooks  such  as  Mohan  et  al.  [37]. 

The  simulation  of  aMOSFET’s  gate-voltage  step  is  illustrated 
in  Fig.  11.  The  MOSFET  was  modelled  according  to  a  data 
sheet  for  a  IRF1404  power  MOSFET  provided  by  International 
Rectifier,  Inc. 

5.5.  Effects  of  parameters  on  the  dynamics 

Examining  Eq.  (24),  it  is  possible  to  understand  the  influence 
of  some  parameters  on  the  dynamic  response  of  a  fuel  cell. 

The  time  constants  are  directly  proportional  to  the  electrode 
capacitance  C.  As  this  capacitance  may  vary  [32],  the  time 
constants  may  vary  accordingly:  the  larger  the  capacitance,  the 
slower  the  transient. 

The  denominator  of  Eq.  (24)  is  made  up  of  two  terms.  When 
both  are  small,  the  transient  will  be  slower;  if,  instead,  one  or  both 
terms  in  the  denominator  are  large,  the  dynamics  will  be  faster. 
This  means  that  the  time  constants  will  be  sensibly  larger  when 
there  is  a  large  sum  of  internal  and  external  resistances  and  a  large 
variation  of  the  overvoltage  with  reaction  current.  Conversely, 
if  any  of  these  two  is  small,  the  time  constants  will  be  smaller. 
This  is  consistent  with  the  observation  of  faster  transients  at  high 
values  of  current,  as  in  Figs.  4  and  5. 

The  time  constants  can  also  be  interpreted  graphically.  In  Fig. 
5  it  was  observed  that  the  electrochemical  transients  that  oc¬ 
curred  when  stepping  the  value  of  the  outer  circuit’s  resistance 
varied  in  duration  depending  on  the  direction  of  the  step.  This 
has  been  confirmed  by  the  model  in  Fig.  10.  It  is  possible  to  visu¬ 
alise  the  cause  of  this  behaviour  by  means  of  the  instantaneous 
characteristic,  or  iso-?;  line,  that  was  previously  introduced. 

The  transient’s  differential  Eq.  (16)  states  that  f  oc  (i  +  ic  — 
?r).  It  is  therefore  important  to  be  able  to  visualise  the  position 
of  i  and  ir  on  the  plot  during  the  transient.  As  long  as  the  load 
is  purely  resistive,  the  operating  point  is  bound  to  be  on  the 
external  circuit  characteristic.  However,  it  must  also  be  on  one 
of  the  feasible  points  for  the  fuel  cell,  which  are  represented  by 


Current,  A  cm'2 


Fig.  1 1 .  The  simulated  transient  resulting  from  having  changed  the  gate  voltage 
of  a  MOSFET  from  3.05  to  3.1  V,  when  this  MOSFET  is  operating  as  the  load 
of  a  PBI  fuel  cell  at  150  °C.  The  markers  are  spaced  by  0.01  s,  except  for  the 
switching  instant  when  there  is  a  spacing  of  1  (xs. 


the  instantaneous  characteristic:  therefore,  the  current  density  i 
always  lies  at  the  intersection  between  the  load’s  characteristic 
and  the  cell’s  instantaneous  characteristic. 

is  the  reaction  current  density  corresponding  to  overvoltage 
ip,  at  any  time,  it  represents  the  rate  of  the  electrochemical  reac¬ 
tion.  The  instantaneous  characteristic  represents  all  the  operating 
points  to  which  the  fuel  cell  can  be  moved  instantly;  its  intersec¬ 
tion  with  the  steady-state  polarisation  curve  represents  therefore 
the  steady-state  point  corresponding  to  the  present  value  of  the 
overvoltage  p,  and  is  at  the  coordinate  i,-  —  ic  on  the  i  axis. 

Therefore,  i  is  found  at  the  intersection  of  the  instantaneous 
characteristic  with  the  external  load,  and  ?,  —  ic  at  the  inter¬ 
section  of  the  instantaneous  characteristic  with  the  polarisation 
curve. 

Looking  at  the  simulated  transient  in  Fig.  10,  it  is  clear  that 
the  drawn  part  of  the  polarisation  curve  is  closer  to  the  char¬ 
acteristic  of  the  higher  resistance,  when  distance  is  measured 
along  the  i  axis.  Since  i]  is  proportional  to  the  distance  between 
the  intersections  of  the  instantaneous  characteristic  with  the  po¬ 
larisation  curve  and  the  external  resistance  along  the  i  axis,  i.e. 
i  +  ic  —  ir,  it  follows  immediately  that  the  transient  to  the  lower 
resistance  will  be  faster  than  the  one  coming  back,  because  of 
the  higher  values  of  i  +  ic  —  iT  along  that  path;  this  is  confirmed 
by  the  experimental  measurements  shown  in  Figs.  4  and  5. 

5.6.  Conditions  for  perfect  power  control 

Considering  the  basic  expression  of  the  voltage  of  a  fuel  cell 
(15),  and  looking  at  Fig.  10,  it  is  evident  that  the  instantaneous 
characteristic  at  a  given  steady-state  point  stays  at  higher  volt¬ 
ages  than  the  polarisation  curve  when  moving  to  higher  values  of 
current  density,  and,  conversely,  at  lower  voltages  for  lower  val¬ 
ues  of  current  density.  This  happens  because  the  instantaneous 
characteristic  is  representative  of  a  single  value  of  which  is 
strictly  increasing  with  current  density.  This  results  in  power 
overshoots  in  any  transient  that  starts  and  ends  in  a  range  of 
currents  lower  than  the  one  associated  to  the  maximum  power 
output,  which  is  often  close  to  the  mass-transfer  limiting  current 


Fig.  12.  A  perfect  response  to  a  step  in  the  set-point  for  power  might  be  obtained 
moving  the  operating  point  to  the  intersection  of  the  iso-??  curve  with  the  iso¬ 
power  curve,  and  following  the  iso-power  curve  until  steady-state  is  achieved. 
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;'l;  it  is  normal  to  assume  that  the  cell  will  operate  below  this 
current  density,  as  operating  above  it  results  in  inefficient  and 
undesirable  operation. 

Under  these  assumptions,  a  perfect  step  change  in  power  out¬ 
put  can  be  obtained  by  leading  the  operating  point  to  the  intersec¬ 
tion  between  the  instantaneous  characteristic  and  the  iso-power 
line  of  the  new  power  output;  this  intersection  always  exists  if 
the  new  power  reference  is  less  than  or  equal  to  the  maximum 
power  that  the  cell  can  produce  at  steady-state.  Then,  by  con¬ 
tinuously  adjusting  the  external  characteristic  (possibly  with  a 
feedback-controlled  MOSFET),  the  operating  point  can  be  led 
along  a  constant-power  line,  until  it  finally  intersects  the  polar¬ 
isation  curve,  where  the  transient  will  stop.  Fig.  12  shows  how 
such  a  transient  should  look  like. 

6.  Conclusion 

A  steady-state  model  of  a  PBI  fuel  cell  was  developed,  and  its 
parameters  were  estimated  by  means  of  parameter  regression  on 
experimental  data.  However,  the  model  fits  well  only  one  branch 
of  the  polarisation  curve,  which  presented  significant  hysteresis. 
The  cause  of  this  hysteresis  is  not  well  understood. 

A  pattern  was  identified  in  transient  tests  of  the  same  cell 
when  changing  the  resistance  of  the  load.  Neglecting  the  slower 
dynamics,  a  dynamic  model  based  on  the  Butler- Volmer  equa¬ 
tion  has  been  developed  for  the  overvoltage  transient.  The  nov¬ 
elty  in  this  model  is  the  use  of  a  function,  the  external  load 
characteristic,  as  input,  rather  than  the  more  common  approach 
of  using  current  or  voltage.  The  model  also  allows  to  use  of 
MOSFETs,  and  any  other  devices  with  varying  characteristic. 
Usage  of  current  or  voltage  as  an  input  would  have  been  un¬ 
fortunate,  as  neither  can  be  a  manipulated  variable  in  a  control 
system. 

This  new  approach  allowed  to  explain  the  observed  transient 
paths  in  a  simple  way,  identifying  the  activation  overvoltage  r]  as 
a  state  of  the  fuel  cell,  together  with  temperature  and  composi¬ 
tions.  The  model  was  also  able  to  simulate  the  different  dynamics 
observed  when  switching  back  and  forth  between  two  external 
resistances.  A  simple  formula  to  estimate  the  time  constants  of 
the  dynamics  has  been  developed,  which  allows  to  predict  an 
approximate  value  for  the  transient’s  time  constant  r  for  an  op¬ 
erating  point. 

Finally,  it  was  shown  how  a  step  in  the  set-point  for  power 
output  can  in  theory  be  perfectly  matched  by  the  power  output 
of  the  fuel  cell,  if  proper  control  could  be  applied,  under  quite 
general  assumptions  on  operating  conditions. 
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